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Abstract 

A highly efRcient, fuhy paralleHzed, six-dimensional tracking model for simulating interactions of colliding 
hadron beams in high energy ring colliders and simulating schemes for mitigating their effects is described. 
The model uses the weak-strong approximation for calculating the head-on interactions when the test beam 
has lower intensity than the other beam, a look-up table for the efficient calculation of long-range beam-beam 
forces, and a self-consistent Poisson solver when both beams have comparable intensities. A performance test 
of the model in a parallel environment is presented. The code is used to calculate beam emittance and beam 
loss in the Tevatron at Fermilab and compared with measurements. We also present results from the studies 
of two schemes proposed to compensate the beam-beam interactions: a) the compensation of long-range 
interactions in the Relativistic Heavy Ion Collider (RHIC) at Brookhaven and the Large Hadron Collider 
(LHC) at CERN with a current-carrying wire, b) the use of a low-energy electron beam to compensate the 
head-on interactions in RHIC. 

Keywords: accelerator physics, parallel computing, beam dynamics 
PACS: 29.27.Bd, 29.27.Fh 



1. Introduction 

In high energy storage-ring colliders, the beam-beam interactions are known to cause emittance growth 
and a reduction of beam life time, and to limit the collider luminosity [1-7]. It has been a key issue in 
a high energy collider to simulate the beam-beam interaction accurately and to mitigate the interaction 
effects. A beam-beam simulation code BBSIM has been developed at Fermilab over the past few years to 
study the effects of the machine nonlinearities and the beam-beam interactions (S-llj. The code is under 



continuous development with the emphasis being on including the important details of an accelerator and 
the ability to reproduce observations in diagnostic devices. At present, the code can be used to calculate 
tune footprints, dynamic apertures, beam transfer functions, frequency diffusion maps, action diffusion 
coefficients, emittance growth, and beam lifetime. Calculation of the last two quantities over the long time 
scales of interest is time consuming even with modern computer technology. In order to run efficiently on 
a multiprocessor system, the resulting model was irnplemented by using parallel libraries which are MPI 
(inter-processor Message Passing Interface standard) Il2ll. state-of-the-art parallel solver libraries (Portable, 
Extensible Toolkit for Scientific Calculation, PETSc) [13|, and HDF5 (Hierarchical Data Format) ^14] . 

The organization of the paper is as follows: The physical model used in the simulation code is described 
in Section [5] The parallelization algorithm and performance are described in Section |31 Some applications 
are presented for the Tevatron, the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider 
(LHC) in Sectional Section [S] summarizes our results. 
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2. Physical model 



In a collider simulation, the two beams moving in opposite direction are represented by macroparticles. 
The macroparticles are generated with the same charge to mass ratio as the particles in the accelerator. 
The number of macroparticles chosen is much less than the bunch intensity of the beam because it becomes 
prohibitive to follow approximately 10^^ particles for even a few revolutions around the accelerator using 
modern supercomputers. These macroparticles are generated and loaded with an initial distribution chosen 
for the specific simulation purpose. As an example, a six-dimensional Gaussian distribution is used for long- 
term beam evolution. The transverse and longitudinal motion of particles is calculated by a sequence of 
linear and nonlinear transfer maps. During the beam transport, a particle is removed from the distribution 
if it reaches a predefined boundary in transverse or longitudinal direction. In our simulation model, the 
following effects are included: head-on and long-range beam-beam interactions, fields of a current-carrying 
wire and an electron lens, multipole errors in quadrupole magnets in interaction regions, sextupoles for 
chromaticity correction, ac dipole, resistive wall wake, tune modulation, noise in lattice elements, single 
and multiple harmonic rf cavities, and crab cavities. The finite bunch length effect of the beam-beam 
interactions is considered by slicing the beam into several chunks in the longitudinal direction and then 
applying a synchro-beam map (l5| . Each slice in a beam interacts with slices in the other beam in turn at 
a collision point. In the following, linear and nonlinear tracking models are described in detail. 

2.1. Transport through an arc 

The six-dimensional coordinates of a test particle in the accelerator's coordinate frame are: x = 
I a;, a; , y, y , z, (5 J , where x and y are horizontal and vertical coordinates, x' and y' the trajectory slopes 



of the coordinates, z = — cAi the longitudinal distance from the synchronous particle, and 5 = /^pz/po 
the relative momentum deviation from the synchronous energy [lf>|. The transverse linear transformation 
between two elements denoted by i and j can be written as 
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(1) 



Here, M. is a. coupled transverse map of off-momentum motion defined by = TZjAii^jTZf^ , where Aii^j 
is the uncoupled linear map described by Twiss functions at i and j elements, and the transverse coupling 
matrix TZ is defined as 

' ( ' (2) 

where is the 2x2 matrix and the symplectic conjugate of the coupling matrix C. The 4x2 dispersion 
matrix is defined hy V — (0,D), and the dispersion vector D = (^Dx, D'^, Dy, Dy^ is characterized by 
the transverse dispersion functions and the map A4, i.e., D = Dj — A^D^ where D^, Dj are the dispersion 
vectors at Since the transport matrix has to be symplectic, the matrix A in Eq. ([T]) is given by 

A — —Ty^S^M.^ where S* is a rearranging matrix (see subsection 12. 7p . The longitudinal map C is given by 

C^[l -(^'/f)^^), where 77 is the slip factor, /3 = w/c, and As the longitudinal distance between the 

two elements, i.e.. As = Sj — Si. It is noted that s is the axis along the beam direction. The nonlinearity of 
synchrotron oscillations is applied by adding the longitudinal momentum change at a rf cavity: 

eVrf 

AS = {sin krfz - sin (l)s) (3) 

where Vrf is the voltage of rf cavity, (j)s the phase angle for a synchronous particle with respect to the rf 
wave, and fc^/ the wave number of the rf cavity. If there are higher harmonic cavities, their effects are added 
to the momentum change. 
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2.2. Beam-beam interactions 

In order to achieve high luminosity in a colUder one can increase the number of bunches which reduces 
the bunch spacing. More bunches can increase the number of parasitic encounters in the interaction regions. 
Since the calculation of beam-beam forces requires large amounts of computational resources, it has to be 
executed rapidly and accurately. BBSIM has three different models for this purpose: a weak-strong model 
for head-on interactions, a look-up table model for long-range interactions, and a Poisson solver model for 
the head-on interactions when both beams have comparable intensities ("strong-strong" model). 

2.2.1. Weak-strong model 

In the weak-strong model we assume that the "weak" beam is affected by the head-on and long-range 
interactions while the opposing beam or "strong" beam is unaffected. The charge distribution of the strong 
beam is assumed to be Gaussian: 

Here, is the number of particles per bunch and q is the charge per particle. Note that the coordinates 
(a;, y, z) are measured in the rest frame of the strong beam. The beam-beam force between two beams with 
transverse Gaussian distribution p (x, y) = f dzp (x, y, z) is well-known [l8| . and the expression for the slope 
change is given by, for elliptical beam with <Jx > <^y'. 

2Nro ( Im [F (x, y)] \ 

where 

F{x,y) = w I ' : I - e ^'1 w I I . (6) 





Here, w (z) is the complex error function defined by w {z) ~ e ^ \ ^ J ' ^^"^ Lorentz 

factor. The constant tq is defined as tq = qq^,/ATTeomoc^, where is the electric charge of the test particle, 
and mo the rest mass of the particle. 

2.2.2. Look-up table model 

The charge distribution of the strong beam in the weak-strong model is not varied during the simulations. 
It is redundant to re-calculate the beam-beam force at every parasitic location and every turn. A look-up 
table is one way to avoid it. The look-up table is used to replace a run time computation with an array 
indexing operation. The beam-beam force of a Gaussian beam distribution is described by the complex 
error function, as shown in Eq. ^ . The calculation of the complex error function can substantially slow the 
beam-beam simulation. However, the look-up table is pre-calculated and stored in a memory, usually in an 
array. When the value of the error function is required, it can be retrieved from the table by an interpolation 
scheme, instead of using Eq. The look-up table method can significantly reduce a computational cost. 
The property of the complex error functions yields the symmetry relations of function F (z) as 



F i-z) = -F (z) , F{z) = -F (z), F (-z) = F (z) (7) 

where z = a: -I- is a complex variable. The symmetry conditions of the function F (z) can reduce memory 
space to store the function values. It is sufficient to build the table for the values of function F (z) in the 
first quadrant of the complex plane, i.e., > and \y\ > 0. 

Interpolation techniques are required to predict a value of a function at a point inside its domain based 
upon the known tabulated values. For a given set of data points (zi, fi), i = 0, . . . , N, where no two z^'s are 
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the same, the interpolated value g (z) at a value z ^ Zi is found from 

N 

= (8) 

i=0 

where the Li is Lagrange's A^-th order polynomials 

N 

u{z)= n f^f- 

In order to save the interpolation time further, one can divide z-space and apply a different degree of the 
Lagrange polynomial. For an example, we apply a sixth order polynomial for small amplitudes \z\ < Aa 
while a third order polynomial is applied for \z\ > 4cr, because the function F (z) varies more rapidly at 
small l^l and slowly at large \z\ . 

2.2.3. Poisson solver model 

The weak-strong model is a good approximation when one beam has much smaller intensity than the 
other, but it is not valid when the intensities of the two beams are comparable because each beam's param- 
eters are changed by the other beam. One has to solve for the field of each beam self-consistently. The fields 
are the solutions of the Poisson equation given by 

V^c?!) (r) = -47rp (r) (10) 

where is the electrostatic potential and p the density function of the beam. The solution can be obtained 

by 

0(r) = J G(r,ri)p(ri)dri (11) 

where G is the Green's function of the Poisson equation and in two space dimension, is given by 

1 

in 

Equation (jlip can be efficiently calculated using a convolution theorem and inverse Fourier transform: 



G{x,y: xi,yi) = -^\n {x - xif + [y - yif . (12) 



0(r)=. J-i(G(a;)p(u;)j (13) 

where G (w) = (^^) ^ J^^ G (r) e-"^ ''dr and p (a;) (^^) ^ J^^ p (r) e'^^ ^dr. It is assumed in Eq. ^ 

that the density function p (r) is periodic in both x and y directions. However, since the beam has a finite 
charge distribution surrounded by a conducting wall in an accelerator system, the transverse beam density 
does not meet the periodicity requirement of EFT techniques. In order to apply the above formalism, the 
density function should be rewritten by, in the doubled computational domain [19|: 

P„e.(x,y) = ^^"'^) ,0<x<L., 0<y<L, ^^^^ 
1^0 , Lx < X < 2Lx, or Ly < y < 2Ly. 

Green's function is defined in the doubled domain, as follows: 

'G{x,y) ,0<x<L^,0<y<Ly 

G„e» (X, y)^{^ ^'"^^ - ,L.<x<2L,^,0<y<Ly 
nen,y ,y, \Q^^ 2Ly-y) , < x < L,, Ly < y < 2Ly ^ ^ 

^G{2Lx - X, 2Ly - y) , L^ < x < 2Lx, Ly < y < 2Ly. 

Both pnew and Gnew are doubly periodic functions with periods 2Lx and 2Ly. It is noted that only the 
potential within a domain {0,Lx] x {0,Ly] is valid. The potential outside the domain is incorrect, but it 
doesn't matter because the physical domain of interest is (0, L^] x (0, Ly]. When one beam is separated far 
from the other, one can apply a shifted Green's function approach [20j. 
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Figure 1: Definition of crossing angles a and (j): a is the crossing plane angle in the x — y plane and 4> is 
the half crossing angle in the x — s plane, s is the axis along the beam direction when there is no crossing 
angle. The x — s plane is the crossing plane defined by the angle a. The beam trajectories, shown by lines 
with arrows, lie in the crossing plane. 



2.2.^. Crossing angle 

When there exists a finite crossing angle between two colliding beams at an interaction point, the beam- 
beam force experienced by a test particle will have transverse and longitudinal components because the 
electric field generated by the opposing beam is not perpendicular to the particle velocity anymore. The 
existence of a longitudinal force makes it difficult to apply the result of previous sections. A transformation 
can be used to remedy the difficulty. It transforms a crossing angle collision in the laboratory frame to a head- 
on collision in the rotated and boosted frame which is called the head-on frame [21I . The transformation 
can be described by a transformation from the accelerator coordinates to Cartesian coordinates, a Lorentz 
boost, and again a backward transformation to the accelerator coordinates: 

X* — z cos a tan + a; [1 + /i* cos a sin 0] + yh*^ sin a sin 
y* = z sin a tan (j) + y [l + h* sin a sin 0] + xh* cos a sin (j) 

z* — + /i* \x cos a sin 6 + y sin a sin 0] 

cos 

^ -px , tan(/) (16) 

jO^ — — /iCOSQ! 



COS <p COS <p 

Py tan0 

p,, — — /ism a 

" cos (p cos (p 

p* — Pz — Px cos a tan (f) — py sin a tan (f> + h tan^ (p 
where a star (*) stands for a dynamical variable in the head-on frame, the Hamiltonian h {px,Py,Pz) 



Pz + l- ^JiPz + l) -pI-pI, K = dh*/dpl, h* {pl,p*y,pl) = h {p*,p*y,pl), a the crossing plane angle in 
the X — y plane, and (f) the half crossing angle in the x — s plane as shown in Fig. [T] 

2.3. Finite bunch length 

The effects due to the finite (as opposed to infinitesimal) bunch length need to be considered when the 
transverse beta functions at the interaction point are small and comparable to cr^. The finite longitudinal 
length is considered by dividing the beam into longitudinal slices and by a so called synchro-beam map flH'l . 
We make slices of both beams moving in opposite directions. Each slice of the strong bunch is integrated 
over its length, and has only a transverse charge distribution at its center. We take into account the collision 
between a pair of slices: the i*^^ slice of a bunch and the j*^ slice of a bunch in the other beam. The collision 
takes place at collision point S (^z^, zi^ = ^ ^z' — zi^ which is usually different from the interaction point. 
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For example, the slice of a bunch has successive collisions with slices of a bunch in the other beam. In 
addition, the electric field varies along the bunch due to the inhomogeneity of the charge density in the 
longitudinal direction, and couples transverse and longitudinal motions. The coupling can be modeled by 
the synchro-beam map which includes beam-beam interactions due to the longitudinal component of the 
electric field as well as the transverse components. The transformation is given by [l^ 



x + S {z, Zf) 
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: Px 
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dy 

= z, 5""'" = 6 - 



dU 

new _ „ _ _ 



1 dU 




1 dU 




1 dU 




1 dU 




1 dU 


2 dx 


s 


~ 2 dx 


s. 


~ 2 dy 


s 


.^^ " 2 dy 


S- 


~ 2 dz 



Here, \g represents the evaluation at the collision point S{z,z 
U = q^/Eq and is given by 



U (x, y; (s) , ay (s)) = 



1 



(17) 

U is the normalized potential energy 
l + exp(-5^ 



y 



^(2a2+C) (2CT2 + C) 



(18) 



The dependence on the bunch length is contained in (7x{s), o'a(s)- The transverse derivatives of the potential 
energy are 

^-Ay'iX,Y;S{z,z,)) (19) 



dU 
dx 



= -Ax'{X,Y;S{z,z^)), — 
s 9y 



where {X,Y) are the transverse coordinates at S'(z,z*), andAx' and Ay' are given by Eq. ([5|). The 
longitudinal derivative of the potential energy which is related to the longitudinal beam-beam kicks is 
expressed by 



dU 
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2N^ro ( ay 
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(20) 



7 V f^y ' 



Note that and have zero amplitude and change their sign at the interaction point if = cty 



0. 



Test particles experience longitudinal acceleration and deceleration passing through the bunch moving in 
the opposite direction. 



2.4- Compensation schemes 

In storage-ring colliders, a beam experiences periodic perturbations when it meets the counter-rotating 
beam in a common beam pipe. The head-on beam-beam interactions occur when the beams collide in the 
detectors while the long-range interactions occur when the beams are simultaneously present at the same 
location but are separated transversely. The nonlinear forces due to these beam-beam interactions result 
in a tune spread and can cause emittance growth, a reduction of beam life time, and therefore reduce the 
collider luminosity. The combination of beam-beam and machine nonlinearities excites betatron resonances 
which can cause particles to diffuse into the tails of the beam distribution and even to the physical aperture. 
Different compensation methods have been proposed: a current-carrying wire for the effects of the long- 
range interactions [22j and an electron lens for the head-on interactions in proton machines ^24-26] . Beam 
collisions with a crossing angle at the interaction point are often necessary in colliders to reduce the effects 
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of the long-range interactions. The crossing angle reduces the geometrical overlap of the beams and hence 
the luminosity. A deflecting mode cavity, also known as a crab cavity, offers a promising way to compensate 
the crossing angle and to realize effective head-on collisions (27l . [28| . We now describe the modeling of these 
compensation schemes in the program. 

2.4- 1- Current- carrying wire 

When the separations at long-range interactions are large compared to the rms beam size the strength 
of these interactions is inversely proportional to the distance. Its effect on a beam can be compensated by a 
current-carrying wire which creates a magnetic field with the same - dependence. This approach is simple 
and it is possible to deal with all multipole orders at once. For a finite length Z^, embedded in the middle of 
a drift length L, the transfer map of a wire can be obtained by 

Ml^^ ^DL/2oMi''^ oDl/2 (21) 

where I?l/2 is the drift map with a length -j, and M^^'' is the wire kick integrated over a drift length. This 
kick map A^^^'' is reproduced by the following changes in slope [13 

f Ax' \ _ Iwlw u-v f X \ , . 

\ Ay' J - 47r(Bp):r2+y2 [y J ^ ' 

where 1^ is the current of the wire , u = \J + IwY + + y'^ and v = \J — Iw)^ + x^ + y^. We also 
take into account the wire misalignment including pitch and yaw angles (O^, 9y) respectively as well as lateral 
shifts (Ax, Ay). The transfer map of a wire can be written as 

Mw = Sax.Av ° ^ij'fiy ° ° M\^^ o Dl/2 ° Tg^fi^ (23) 

where Tg^^g^ represents the tilt of the coordinate system by horizontal and vertical angles 9x,dy to orient 
the coordinate system parallel to the wire, and SAx.Ay represents a shift of the coordinate axes to make 
the coordinate systems after and before the wire agree. When the wire is parallel to the beam, Eq. (|23p 
becomes Eq. (j2ip . For canceling the long-range beam-beam interactions of the round beam with the wire, 
one can get the desired wire current and length by equating Eq. ((22|) and Eq. ([5]); the integrated strength 
of the wire compensator is related to the integrated current of the beam bunch as I^lw — cqN. 



2.4-2. Electron lens 

For the head-on proton-proton beam collisions, particles of one proton bunch are focused by a space 
charge of the counter-rotating proton bunch. The beam-beam effect on the particles of the proton bunch 
can be compensated by a counter-rotating beam of negatively charged particles, for example, a low-energy 
electron beam. In order to cancel out the transverse kick by the counter-rotating proton bunch, the electron 
beam should have the same transverse charge profile and current as the proton bunch. The proton bunch 
typically exhibits an approximately Gaussian transverse profile. If we choose a Gaussian distribution of the 
electron beam, the transverse kick on particles of the proton bunch from the electron beam is given by 



Ax' 
Ay' 



—^C{x,y.(Te) 



(24) 



where N^. is the number of electrons of the electron beam adjusted by the electron beam speed, tq the classic 



proton radius, 7 the Lorentz factor, r^ 
The function ^ is given by 



y 



2, and (7„ the transverse beam size of the electron beam. 



C,{x,y ■■ (Te) 



1 — exp 



2ct.. 



(25) 
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For a non-Gaussian electron charge distribution we implement a flat top profile with smooth edges that gen 
erates a linear beam-beam force near the beam center. This flat top beam profile pe (r) = po/ ^1 + (^/'^e 
delivers the transverse kicks given by Eq. but the function C is as follows: 



c 



1 



-I- tan ^ 9+ + tan ^ ( 



(26) 



where p is a constant, and 9± — \f2 ( — ) ± 1- 



24.3. Crab cavity 

When a particle passes through a crab cavity structure, it experiences a transverse deflection and a small 
change in its longitudinal energy. Crab cavities can compensate for the horizontal or vertical crossing angle 
at the interaction point by delivering oppositely directed transverse kicks to the head and the tail of the 
bunches. In the case of a horizontal crossing, the kicks from the crab cavity are given by 



Ax = -— sm 0s H , Ad = -— cos {(ps ~\ • - 

i^o^ c/ ho \ c / c 



qV 



(27) 



where q denotes the particle charge, V the voltage of crab cavity, Eq the particle energy, (ps the phase of 
the synchronous particle with respect to the crab-cavity rf wave, uj the angular frequency of the crab cavity, 
c the speed of light, z the longitudinal coordinate of the particle with respect to the bunch center, and x 
the horizontal coordinate. In general this is a nonlinear map which introduces synchro-betatron coupling, 
but for small z, this reduces to a linear map in the horizontal-longitudinal plane. The crab cavity causes a 
closed orbit distortion dependent on the longitudinal position of particles, and the beam envelope is tilted 
all around the ring. For a bunch shorter than the rf wavelength of the crab cavity deflecting mode, the tilt 
angle of the beam envelope at a location with a beam position monitor (BPM) is given by 



tan 0c 



crab 



c^po 



cos {Aip — ttQ) 



2 sin ttQ 



(28) 



where /? is the beta function at the BPM position, /3crab the beta function at the crab cavity, Acp the phase 
advance between the crab cavity location and the BPM, and Q the betatron tune. The simulations of a 
crab cavity in the SPS accelerator at CERN using BBSIM will be described in another paper. 



2.5. Particle distribution 

At the beginning of a simulation, the simulation particles are distributed over the phase space x = 
{x,x' ,y,y' , z,6)'^ , called the initial loading. In any simulation the number of particles N is limited by the 
computational power. In order to make the best use of a small number of simulation particles compared to 
the real number of particles in the accelerator, the loading should be optimized. Indeed the initial loading 
is very important because this choice can reduce the statistical noise in the physical quantities. 

Gaussian distribution: For long-term particle tracking where we calculate emittance growth, we consider 
an exponential distribution in action (Gaussian distribution in coordinates) of the form: 



P (x) = po cxp - 



J X 



J 11 



J. 
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where J^, Jj,, and Jz are the transverse and longitudinal action variables deflned by 



(29) 



x"^ + [(3xX + axX 



1 



8 Rv, 
TT \r]\ 



(30) 



where R is the radius of the accelerator, h the harmonic number, i/g the longitudinal tune, E and K the 
complete elliptical integrals, and 

,2^2 / 2 



CTj^ , a J , and ctj^ are the rms sizes of action variables. The simulation particles are generated by two 
steps: 

1. The action variables {Jx, Jy, Jz) of particles can be directly generated from the distribution function 
by the inverse transform method and the bit-reversed sequence [30]. 

2. For example, x and x' are correlated and their distribution is p{x,x') — po exp ^— ^ +{/3:^x^+a^x) 

Since the horizontal action is determined at the first step, the horizontal coordinates {x,x') can be 
obtained from the random variates: 

X = \/Jx cosdx, x' = ^/jx {ainOx ~ ax cos9x) / Px 
where the value of Ox is randomly distributed within the interval < Ox < 27r. 

Hollow Gaussian distribution: In most cases of particle tracking, lost particles are observed only above 
a certain large transverse action while the beam core is stable. An example is shown in Section 14.11 A 
hollow beam is a beam with zero central intensity along the longitudinal beam axis. For the generation of a 
hollow beam, a bunched beam distribution in longitudinal phase space is a Gaussian, but a distribution in 
transverse phase space is a hollow Gaussian. The procedure of generating the hollow distribution is the same 
as that for the Gaussian distribution except that the amplitude of transverse action of a particle should be 
larger than a minimum value, i.e., Jx + Jy ^ crj. Since most of the stable particles are not included in the 
tracking simulation, the hollow beam model simulates a large transverse amplitude Gaussian distribution 
using a small number of macro-particles. This distribution is useful when calculating beam lifetimes. 

2.6. Particle diffusion 

Diffusion coefficients can characterize the effects of the nonlinearities present in an accelerator, and can 
be used to find numerical solutions of a diffusion equation for the density 31, 3^. The solutions yield the 



time evolution of the beam density distribution function for a given set of machine and beam parameters. 
This technique enables us to follow the beam intensity and emittance growth for the duration of a luminosity 
store, something that is not feasible with direct particle tracking. The transverse diffusion coefficients can 
be calculated numerically from 

A, (a^,a,) = ^ {{J^{a,,N) - J,(a„0)) iJj{a,,N) - Jj{a,,0))) (32) 

where Ji (a^, 0) is the initial action at an amplitude a^, Ji (a^, N) the action with initial amplitude after N 
turns, the average over simulation particles, and (i, j) are the horizontal x or the vertical y coordinates. 
Equation (j29p is averaged over a certain number of turns to eliminate the fluctuation in action due to the 
phase space structure, e.g. resonance islands. These diffusion coefficients can be directly used to compare 
amplitude growth under different circumstances, e.g with different tunes. Emittance growth and beam 
lifetimes can be calculated when these coefficients are used in a diffusion equation, as mentioned above. 

2. 7. Symplecticity 

In the absence of dissipative effects, particle motion in an accelerator can be described by Hamilton's 
equations of motion. Hamiltonian systems obey the symplectic condition which guarantees the conservation 
of phase space volume as the system evolves, this is also known as Liouville's theorem. For transfer maps 
described in previous subsections the symplectic condition requires 



M^SM = S, 5=0 s (33) 




where s = 



( 



1 
-1 



) 



is an antisymmetric 2x2 matrix, and M is a transfer matrix for a linear system or 



the Jacobian matrix of a nonhnear map around any particle trajectory. For a nonlinear map M : x — > x, 
the Jacobian matrix is obtained from first-order partial derivatives of the new coordinates with respect to 
the old ones. The elements are defined as My — dxi/dxj. During implementation of the maps for beam 
dynamics, one should check to ensure that the map is as symplectic as possible. As a measure of the 
symplecticity, a matrix norm of || M-^S'M — "S*]! is used in BBSIM. The accuracy of the look-up table model 
mentioned in subsection I2.2.2[ for example, depends on the number of sample points in a given complex 
space needed for interpolating the function. Poor interpolation accuracy may violate the symplecticity. and 
lead to emittance blow-up or shrinkage. We use the symplectic norm obtained with the direct calculation of 
the complex error function as the benchmark. We find for example, that in order to maintain the symplectic 
norm with 70 long-range beam-beam interactions in the Tevatron, the number of sample points should be 
more than 4 points per rms beam size. 

2.8. Diagnostics 

Numerical simulation enables the generation of very large amounts of data. The BBSIM code monitors 
physical quantities, for example, particle amplitudes and saves them into an external file during the sim- 
ulation. According to a problem of interest, the quantities to be saved can be chosen in order to extract 
valuable information from post-processing. In addition, some diagnostic functions are calculated in the code 
as follows: 

Betatron tune distribution: The betatron tune in an accelerator is one of the most important beam 
parameters. The tune of each particle in the beam distribution is calculated with a Hanning filter applied 
to an fast- Fourier transform of particle coordinates found from tracking [SS^ . 

Beam transfer function: The beam transfer function (BTF) is defined as the beam response to a small 
external longitudinal or transverse excitation at a given frequency. BTF diagnostics are widely employed in 
accelerators due to its non-destructive nature. A stripline kicker or rf cavity excites betatron or synchrotron 
oscillations respectively over the appropriate tune spectrum. The beam response is observed in a downstream 
pickup. The fundamental applications of BTF are to measure the transverse tune and tune distribution 
by exciting betatron oscillation, to analyze the beam stability limits, and to determine the impedance 
characteristics of the chamber wall, and feedback system |34i] . In the code, we apply a sinusoidal driving 
force to a beam in a transverse plane. The driving frequency is swept in equidistant steps over a continuous 
frequency range which includes the betatron tune. At each new frequency there is initially a transient 
response which must be allowed to relax before the frequency is extracted from the data. We avoid the issue 
of the transients in the simulations by reloading the initial particle distribution at each new frequency. 

Frequency diffusion: We have calculated frequency diffusion maps as another way to investigate the 
effects of nonlinear forces. The map represents the variation of the betatron tunes over two successive sets 



of the tunes [35]: The variation can be quantified by c? = log \/^'^x + ^^y' where {Aiy^ = Vx —Vx \ l^Vy — 



Vy — Vy ) are the tune variations between the first set and next set of 1024 turns. If the tunes [v^ ,v\ 



are different from \ Vx \ vy j, the particle is moving to different amplitudes. A large tune variation is 
generally an indicator of fast diffusion and reduced stability. 

Dynamic aperture: The dynamic aperture of an accelerator is defined as the smallest radial amplitude of 
particles that survive up to a certain time interval, for example, 10^ turns. As the number of turns increases, 
the dynamic aperture approaches an asymptotic value. Initial particles are distributed uniformly over the 
transverse phase space with amplitudes typically varying between 0-20 cr, where a is the rms transverse 
beam size. The longitudinal amplitude is chosen as largest value within a bunch. 

Emittance: The emittance is defined as the area (or volume) of phase space enclosed by the ellipse 
containing all the particles in its interior. Statistically, the rms beam emittance can be calculated by a 
determinant of S-matrix of a beam distribution: 






[det(E)] 



(34) 
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where d is the dimension of phase space, the element of S-matrix is — {{d — (Q)) {Cj ~ (Cj)))i ^-'id 



addition to the emittance of each degree of freedom, four- and six-dimensional emittances are calculated to 
see the correlation and coupling between the phase space coordinates. 

Beam loss: The beam loss is one of the fundamental observables and it can be directly compared with 
simulation. During a beam simulation, each particle is monitored if it reaches a physical aperture transversely 
or the rf bucket longitudinally. The particle passing beyond the aperture is considered as a lost particle. 
Unlike a real machine, several virtual apertures are placed inside a beam pipe. The multiple apertures are 
used to find beam losses at different apertures. 

3. Parallelization 

Realistic simulations of beam dynamics demand large computational resources. Calculations on these 
large number of particles can be distributed over several processors of a parallel computer to improve 
performance. Two basic approaches exist to allocate the calculations to the processors, particle based and 
domain (space) based partitions. In the former approach, the particles are uniformly allocated to the 
processors. They are not limited to a certain spatial domain. The completion time of a parallel solution 
depends on the processor with the maximum computational workload. The particle decomposition can 
distribute the computational load evenly among all processors while the interaction between particles, for 
example, intra-beam scattering needs a very large number of communications between processors since 
the interacting particles can be located in a distant processor. Conversely, in the domain decomposition 
approach, the spatial domain is partitioned into elementary regions, and each processor is responsible for one 
of these regions. The particles in the accelerator simulation are transported by the lattice map. The map 
causes significant particle movement which may cause the load to become quickly unbalanced. The simulation 
of colliding beams has two aspects, i.e., pure particle transport and electromagnetic field evaluation. The 
domain deposition approach is an efficient way of parallelizing the field solver. To achieve the workload 
balanced, our approach is to use both decomposition schemes. 

We have implemented a parallel calculation in the BBSIM code to perform a tracking simulation of large 
numbers of particles. When the weak-strong beam-beam model is used, only the particle decomposition 
scheme can be applied for parallel computation. Its implementation can be made trivially because the 
macroparticles are never moved from one processor to another. No inter-processor communication is nec- 
essary while the particle trajectories are being developed. Most calculations on each node are executed 
sequentially. In this model the communication between the parallel processes is only required for reading 
input data, generating an initial beam distribution, calculating diagnostics such as beam emittance, and 
writing out the diagnostic information. For the Poisson solver model, however, we have used a particle-in-cell 
(PIC) model to update the electromagnetic field. The PIC model represents the beam as a large number of 
computational particles moving according to classical mechanics. The PIC algorithm can be characterized 
as follows: (a) integrate over particles to obtain a charge distribution on the grid point, (b) solve a Poisson 
equation for the potential, and (c) interpolate the potential or field onto particles for a small interval of 
time to advance the position and velocity of particles. Part (a) requires O (Ng) numeric operations for a 
EFT Poisson solver, where Ng is the number of grid points per dimension and d is the number of degrees 
of freedom. Part (a) and (c) obviously require O (Np) operations, where Np is the number of computation 
particles. In general, Np is much larger than Ng in that the number of particles should increase according 
to the degree of freedom to maintain the statistical noise to be constant in a higher spatial dimension. The 
particle calculations thus dominate the overall computational process, which suggests a prior parallelization 
of particle calculation. Master/slave configuration of computational nodes shown in Fig. [2] is considered due 
to the difference of numeric operations between particles and field updates. 

Each processor on the master and slave nodes possesses the same number of particles. All processors 
are responsible for advancing their particles. On the contrary, the master node may be a single or many 
processor(s), depending on the number of grid points required. The charge density of a beam is deposited 



= {x, x' , y, y' , z, S}. For example, horizontal emittance is obtained by — det 
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Figure 2: Master/slave communication diagram. 
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Figure 3: Plots of (b) parallel speedup versus the number of nodes, and (b) CPU time versus the number of 
simulation particles, cerf and table represent the weak-strong model, and look-up table model respectively. 



on the computational grids of each processor using standard area weighting (or higher order) methods [36|. 
The master node gathers the charge density from all processors, and solves the Poisson equations in parallel. 
The master node broadcasts the solution of the electric field to all processors such that each processor exerts 
the electromagnetic force on the particles owned by the processor. 

The performance of the master/slave parallelization approach has been investigated using a real lattice 
of the Tevatron which has two head-on beam-beam collisions and 70 long-range beam-beam interactions. 
Speedup test has been performed on the Cray XT5 of the National Energy Research Scientific Computing 
Center at Lawrence Berkeley National Laboratory. The system is built up of 664 nodes with two quad-core 
AMD 2.4 GHz processors per node. The speedup of a parallel program is a measure of the utilization of 
parallel resources and is simply defined as the ratio between sequential execution time and parallel execution 
time : 

5p = ^ (35) 

where -p is the number of processors, T\ is the execution time of the sequential algorithm, and Tp is the 
execution time of the parallel algorithm with p processors. For a fixed number of processors p, typically the 
speedup is < S'p < p. Ideally all parallel programs should exhibit a linear speedup, i.e., Sp = p, but it is 
not common because communication between processors is considerably slower than computation in each 
processor. Figure [2] (a) illustrates the resulting speedup as a function of the number of processors. 

The parallelization speedup based on the total simulation time is compared for simulations with the 
weak-strong model and the look-up table model. The speedup curves are very close to the ideal one below 
a certain number of processors, while they are less than optimal when the number of processors increases 
above a critical value, for example, 2^ processors. On large numbers of processors a relative fraction of the 
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communication time in the total computing time becomes large. A parallel efficiency, defined as the speedup 
factor divided by the number of processors, can be obtained as high as 87% up to the critical number of 
processors. Though the efficiency falls well below 38% when the number of processors is beyond 2^", it runs 
367 times faster than on a single processor. In order to see the scalability of our parallel code for larger 
problem sizes. Fig. [3] (b) shows the execution time as a function of the number of macro-particles. Here the 
number of processors is fixed at 2^ for all cases. It is seen that with increasing the number of simulation 
particles, the execution time also increases linearly. 

4. Applications 

In high energy storage-ring colliders, the beam-beam interactions cause emittance growth, may reduce 
beam lifetime, and hence limit the collider luminosity. We have used BBSIM to study beam-beam interactions 
and their compensations in the Tevatron, in RHIC and in the LHC. 

4..I. Tevatron 

The luminosity of a collider is found from 

C ^ ^f^f^-R (36) 

where A^i and N2 are the bunch populations of the colliding beams, / the revolution frequency, Nb the 
number of bunches in one beam, and ay the horizontal and vertical rms beam sizes at the collision points 
respectively, and R the luminosity reduction factor due to the "hour-glass" effect and due to non-zero crossing 
angle at the interaction point. The beam-beam tune shift of beam 1 is proportional to the factor N2/<Jx'^y 
and experience from colliders worldwide has shown that the achievable tune shift (and hence luminosity) is 
limited by the dynamics of the beam-beam interaction. In the Tevatron, proton and anti-proton bunches 
collide at two detectors called CDF and DO. They share the same beam pipe. Since the two beams circulate 
on helical orbits, the optics and dynamics of the beam-beam interactions are complex. The beam-beam 
interactions occur all around the ring and at varying betatron phases. In run II, each beam has three trains 
of 12 bunches [38]. Each bunch experiences 72 interactions: two interactions are the head-on collisions in 
the detectors. However, the other 70 interactions are long-range, and are placed at different locations for 
each bunch. Consequently the beam separation distances between proton and anti-proton beams at the 
long-range locations are different from bunch to bunch. Figure 3] shows the radial beam separation of three 
anti-proton bunches from the proton bunches in units of the rms beam size of the proton beam at the 
locations of the beam-beam interactions. The long-range interactions of special importance are those on 
either side of the head-on interaction points. These occur at small separations and the beta functions there 
are large. It was observed that the emittance growth at the end bunches of each train is smaller than those 
in the middle of the train. Here we choose two end bunches {^1 and #12) and one middle bunch (#6) of 
the first train. 

Beam emittance growth and loss rate are routinely measured during the Tevatron operation. They can 
be directly compared with numerical simulations but only for relatively short times. Figure [S] (a) shows the 
time evolution of the four-dimensional emittance of bunches #1, #6, and #12 for 15 hours of high energy 
physics (HEP) run of store # 7650. The emittance is calculated and plotted by €4^ = ^t^^y It is observed 
that during the HEP run, the emittance growth is nearly linear. The growth rate is 6.7%/hr. Figure [S](b) 
shows the measured beam loss rates of anti-proton bunches during the first 1 hour of store #7601-#7650 at 
collision energy 960 GeV. In order to see the effects of beam-beam interactions on the beam loss, the loss 
rate is obtained by subtracting the particle losses due to luminosity at the main interaction points from the 
total beam loss rate. Averaged loss rates of bunch #1 and #12 are 1.4 %/hr and 1.2 %/hr respectively, 
while the loss rate of bunch #6 is 2.3 %/hr. We performed the simulations of emittance growth and particle 
loss of anti-proton beam, as shown in Fig. [S] (c)-(d). The particle tracking is carried out over 10 turns 
corresponding to approximately 3.5 minutes storage time of the Tevatron. In the simulation, nominal tune is 
(20.571, 20.569). Initial transverse (95% normalized) emittance of anti-protons (ea;,ej,) is set to be (9.0,7.8) 



13 




16 
14 
12 
10 
8 
6 
4 
2 



T 1 1 1 1 r 



#06 



J I I U I L 



1 2 3 4 5 6 
s [km] 



16 
14 
12 
10 
8 
6 
4 
2 




#12 



1 2 3 4 5 6 7 
s [km] 



Figure 4: Separation distance between proton and anti-proton beams for anti-proton bunches #1, #6 and 
#12. The separation is normahzed by proton beam's rms size. 



mm-mrad from averaging the measured emittances while proton's initial emittance is (18,23) mm-mrad. 
Bunch intensities of anti-proton and proton are 0.86 x 10^^ and 2.64 x 10^^ respectively. Figure[n](c) shows 
the emittance growth of three bunches during the simulation. The growth rate is approximately 9 %/hr, 
which is close to the measured growth rate 7 %/hr in Fig. [S](a). The emittance does not vary from bunch to 
bunch. However, the beam losses vary considerably from bunch to bunch. As shown in Fig. [5](d), bunch #6 
loses more particles than bunches #1 and #12, which agrees well with the observation. For the simulation 
of beam loss, we used the hollow Gaussian distribution in transverse action coordinates. Most of the lost 
particles have large transverse actions as shown in Fig. IH] (a) , while the lost particles are distributed over 
the entire range of longitudinal action, as shown in Fig. [S] (b). 

The compensation of long-range effects in the Tevatron with a current-carrying wire was investigated 
using an earlier version of the code [s^ . It was found that a single wire was unable to compensate for all the 
70 interactions, since they were all at different betatron phases from the wire. 

4-2. Relativistic Heavy Ion Collider 

We have studied the effects of a current-carrying wire on the beam dynamics in RHIC [3^ . Two current- 
carrying wires, one for each beam, have been installed between the magnets Q3 and Q4 of IP6 in the RHIC 
tunnel. In the physics run 9, an attempt was made to compensate the long range beam-beam interaction 
which shows the reduction of beam loss [s^. During the physics run 7 and 8, the impact of current-carrying 
wires on a beam was measured without an attempt to compensate the beam-beam interactions. However, 
the experimental results help to understand the beam-beam effects because the wire force is similar to the 
long-range beam-beam force at large separations. As an example. Fig. [7] plots the beam loss rate due to the 
wire as a function of beam-wire separation distance. The onset of beam losses is observed at 8 cr and 9 a 
for gold and deuteron beams, respectively. The threshold separation for the onset of sharp losses observed 
in the measurements and simulations agree to better than 1 cr. It is also significant that the simulated loss 
rates at 7 and 8 a separation for the gold beam and 8 and 9 cr for the deuteron beam are very close to the 
measured loss rates. At fixed separation, the wire causes a much higher beam loss with the deuteron beam 
than with the gold beam. The loss rate for the gold beam at a 8 cr separation is about 10 %/hr while for the 
deuteron beam the loss rate is about an order of magnitude higher both in measurements and simulation. 
Simulations of the beam loss rate when the wire is present are in good agreement with the experimental 
observations. 

In the proton-proton runs of RHIC, the maximum beam-beam parameter reached so far is about — 
0.008. This tune shift is large enough that the combination of beam-beam and machine nonlinearities excite 
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Figure 5: (a) Variation of anti-proton emittance of three bunches, ^1, ^6, and #12, of store #7650, (b) 
non-luminous loss rates of anti-proton during the first 1 hour of stores #7601-#7650, (c) simulation of 
anti-proton emittance growth, and (d) simulation of anti-proton beam loss. Here the emittance is plotted 
as e4rf = ^e^ey In the simulation, initial anti-proton emittance {ex,^y) is (9.0,7.8) mm-mrad, bunch length 
1.5 usee, and bunch intensity 0.86 x 10^^. Proton's initial emittance is (18,23) mm-mrad, bunch length 1.7 
nsec, bunch intensity 2.64 x 10^^. Nominal tune is (20.571, 20.569). Revolution frequency is 47.7 kHz. 
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Figure 6: (a) Scatter plot of lost particles in action space {yfJZ, ^fJy) '"^"^^ (b) plot of lost particles versus 
\J Jx + Jy for different longitudinal action. The axis variables are normalized by rms size of transverse 
action. 




Figure 7: Comparison of the simulated beam loss rates with the measured as a function of separations, (a) 
gold beam at collision energy, (b) deuteron beam at collision energy ^] . 
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Table 1: Comparison of particle loss for different electron beam profiles and intensities. 



betatron resonances which cause emittance growth and diffuse particles into the tail of beam distribution 
and beyond. Consequently RHIC is actively developing an electron lens for compensating the head-on 
interactions [40]. In order to seek the electron lens parameters at which the beam life time is improved, we 
choose three different electron beam distribution functions: (a) Ict Gaussian distribution with the same rms 
beam size as that of the proton beam tr, (b) 2(7 Gaussian distribution with rms size twice that of the proton 
beam, and (c) Smooth-edge-flat-top (SEFT) distribution with an edge around at 4 a. When the electron 
beam profile matches the proton beam, the full compression of the tune spread requires the electron beam 
intensity A^e 4 x 10^^ which is defined as the electron beam intensity required for full compensation. Table 
[T] shows the results of particle loss for different intensities with the three electron beam profiles. 

At an intensity A^e = 4 x 10^^, the particle loss is nearly six times the loss without beam-beam compensa- 
tion. The beam lifetime at A^e = 2 x 10^^ however is comparable with that of no beam-beam compensation. 
As the electron beam intensity is decreased, the particle loss decreases significantly, and is reduced to 30% of 
that without beam-beam compensation at N,^ = 0.5 x 10^^. For the 2a Gaussian and SEFT electron beam 
profiles, we calculated particle loss for different electron beam intensities. The upper limits of the electron 
beam intensity for these two distributions are chosen so that peak of the electron profile matches that of 
the full compensation at la Gaussian. For the intensities 2 x 10^"'^ and 4 x 10^^ of 2a Gaussian profile, 
there is a significant reduction in beam loss, for example, below 10% of the particle loss without beam-beam 
compensation when the electron beam intensity is 2 x 10^^. A significant improvement of beam lifetime with 
the SEFT profile is also observed below 8 x 10^^. There is a threshold electron beam intensity below which 
beam life time is increased: 2 x 10^^ for the la Gaussian, 8 x 10^^ for the 2a Gaussian, and 16 x lO^^for the 
SEFT profile. Particle loss is relatively insensitive to electron lens current variations below the threshold 
current with the 2a Gaussian and SEFT profiles. This looser tolerance on the allowed variations in electron 
intensity will allow greater intensity fluctuations and is likely to be beneficial during experiments. 

J^.S. Large Hadron Collider 

As mentioned above, long-range beam-beam interactions cause emittance growth or beam loss in the 
Tevatron and are expected to deteriorate beam quality in the LHC. Increasing the crossing angle to reduce 
their effects has several undesirable effects, the most important of which is a lower luminosity due to the 
smaller geometric overlap. For the LHC, a wire compensation scheme has been proposed to compensate the 
long-range interactions [23] . However, several issues need to be resolved for efficient compensation. With 
the design bunch spacing, there are about 30 long-range interactions on both sides of an interaction point 
(IP). The beam-beam separation distance varies from 6.3 a to 12.6 a. The resulting beam-beam force is 
not identical to that generated by a single or multiple wire(s) but can be closely approximated by the wires. 
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Figure 8: Plot of (a) beam-beam separation at IP 1 and 5 and (b) particle loss according to wire separation 
distance with wire strength 82.8 Am. 



Unlike the Tevatron. the long-range forces in the LHC are all at nearly the same betatron phase and this 
makes the compensation scheme feasible. The wire-beam separation distance is one of the parameters which 
determine the performance of a wire compensator. Figure [5] (a) shows the beam-beam separation distance 
normalized by the transverse rms bunch size. Two counter-rotating beams collide at a vertical crossing angle 
near IPl while they collide at a horizontal crossing angle near IPS. The separations are asymmetric with 
respect to the interaction points. The reference wire-beam separation (9 a) is chosen as the average of beam- 
beam separations. Figure |5] (b) shows the results of particle loss for different wire-beam separations. The 
particle loss saturates at large separation while there is a sharp increase of particle loss at small separation. 
We directly see the minimum particle loss between 0.9 and 1.0 of the reference separation. It reveals that the 
average of beam-beam separations is close to an optimal separation between the wire and the high energy 
bunch. 



5. Summary 

In this paper, an efficient parallel beam simulation model for circular colliders is presented in order to 
study the effects of beam-beam interactions and machine nonlinearities, and the effectiveness of beam-beam 
compensation schemes. We have included the major nonlinearities present in accelerators in our program as 
well as models for several methods to compensate the effects of beam-beam interactions. A particle-domain 
decomposition scheme is implemented with the master/slave configuration to achieve a balanced workload 
in a parallel environment. A performance test of beam-beam interactions indicates that the parallelization 
scheme scales linearly in both the number of processors and the number of particles in the beam. We 
have used the program to study the emittance growth and beam loss of different bunches due to the beam- 
beam interactions in the Tevatron, the compensation of head-on beam-beam interactions with a low energy 
electron beam in RHIC, and the long-range beam-beam compensation using a current-carrying wire in the 
Tevatron, RHIC and the LHC. The pattern of beam losses observed in the Tevatron is reproduced in the 
simulations. In RHIC, simulations of the beam loss rate when the wire is present are in good agreement 
with the experimental observations. We have several predictions from the results of head-on compensation 
in RHIC. For example we find that proton beam life time is increased if the electron beam intensity is kept 
below a threshold intensity. An electron beam wider than the proton beam at the electron lens location 
is found to increase beam life time. The results of LHC simulation with the current carrying wire show 
that the particle loss is minimized when the beam-wire separation is close to the average of beam-beam 
separations. 
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